### This script will produce a BRT model of underlog temperature
### Parameters tc, lr, tf, nt, i and output.dir are defined by a subumission script

# Obtain command line arguments from .sh file

args=(commandArgs(TRUE))

# Evaluate the arguments for use in this script

for(i in 1:length(args)) 
	
	{
		
	eval(parse(text=args[[i]]))
 
	}
		
# Load Elith Source Code

setwd('/home1/99/jc152199/brt/')
source('brt.functions.R.cjsedit.r')

# Load gbm library

library('gbm')

#### Reset work directory

setwd(output.dir)

# Read in data to model

model.data = read.csv('/home1/99/jc152199/underlogdownscale/ModelDataforULBRT.csv',header=T)

#### Remove site-days with Empirical Temperatures above 30 degrees

model.data = model.data[which(model.data$ulmax<35),]

#### Calculate temperature range for empirical air temps

model.data$BRTairRange = model.data$BRTairmax - model.data$BRTairmin

# Randomly shuffle model.data before subsetting to training/testing sets

model.data = model.data[c(sample(c(1:nrow(model.data)), nrow(model.data),replace=FALSE)),]

#### Remove sites with both a log density and circumference

model.data = model.data[which(is.na(model.data$logden)==F & is.na(model.data$logcirc)==F),]

# Randomly sample half of data set to produce a training and a testing set

sample.vector = c(sample(c(1:nrow(model.data)), round(nrow(model.data)*as.numeric(tf),0),replace=FALSE))

### Create training and testing sets

train.data = model.data[c(sample.vector),]
test.data = model.data[-c(sample.vector),]

# Run BRT model using gbm.step, parameters are defined as arguments, indepedent variables are limited to BRT max, min, and range
		
brt.gbm.step = gbm.step(data=train.data,gbm.x = c(14,16:20),gbm.y = 3,max.trees=as.numeric(nt),family = "gaussian",tree.complexity = as.numeric(tc), learning.rate = as.numeric(lr) ,bag.fraction = 0.5, n.folds=10)

###  Save model results

save(brt.gbm.step, file=paste(output.dir,'ULModel.Rdata',sep=''))

# Produce partial dependence plots

gbm.plot(brt.gbm.step,variable.no = 0, nt = brt.gbm.step$n.trees, plot.layout = c(2,3))

# Append preds to model data

train.data$ulmaxpreds = brt.gbm.step$fitted

### Write training data (with preds) and testing data (with preds) to the output directory

write.csv(train.data, file=paste(output.dir,'UL_BRT_Model_Data_Plus_Preds.csv',sep=''),row.names=F)

### Now produce some plots from the test and training data and summarize model outputs

# First Predict onto test data

gbm.predicted = gbm.predict.grids(brt.gbm.step, test.data, want.grids = F, sp.name = "preds")

# Calculate total deviance from test data

total.deviance = sum((test.data[,3]-gbm.predicted)*(test.data[,3]-gbm.predicted))/length(gbm.predicted) # Calculate mean deviance (weighted by number of observations)

### Calculate test data residuals

residuals = test.data[,3]-gbm.predicted

# Make some plots from testing data

### Start .png driver

png(paste(output.dir,'Test_Data_Obs_Versus_Preds.png',sep=''),units='cm',height=12, width=12, res=1000)

#### Plotting limits

lims = range(c(test.data[,3],gbm.predicted),na.rm=T)

#### Plot up testing obs versus preds

plot(test.data[,3],gbm.predicted, xlab = 'Obs Test Data', ylab = 'Test Fitted Values', main = paste('Obs vs Preds - Deviance - ',round(total.deviance,4),sep=''), xlim=lims, ylim=lims, col='olivedrab1')

#### 1:1 line

abline(a=0,b=1, lty=2, lwd=1, col='red')

### LM

lm1 = lm(gbm.predicted~test.data[,3]) # Perform a linear model using the Empirical Data as y-values and the preds (AWAP or microCLIM) as x-values

#### Plot lm relationship

abline(lm1,col='lightseagreen',lty=2, lwd=1)

#### Legend

legend('topleft',legend = c(paste('Adj. r^2 ',substr(summary.lm(lm1)[9],1,4),sep=''),paste('Slope ',round(lm1$coefficients[2],2),sep=''),paste('Intercept ',round(lm1$coefficients[1],2),sep='')),text.col=c('lightseagreen','lightseagreen','lightseagreen'), bty='n')

#### Shut device

dev.off()

### Plot observed test data against residuals

### Open .png driver

png(paste(output.dir,'Test_Data_Obs_Versus_Resids.png',sep=''),units='cm',height=12, width=12, res=1000)

#### Plotting limits

xlims = range(test.data[,3],na.rm=T)
ylims = range(residuals,na.rm=T)

### Plot data

plot(test.data[,3],residuals, xlab = 'Obs. Testing Data', ylab = 'Residual (Obs-Fit)', main = paste('Obs vs Resid - Deviance - ',round(total.deviance,4),sep=''),xlim=xlims, ylim=ylims, col='olivedrab1')

### LM

lm1 = lm(residuals~test.data[,3]) # Perform a linear model using the Empirical Data as y-values and the preds (AWAP or microCLIM) as x-values

#### Plot lm relationship

abline(lm1,col='lightseagreen',lty=2, lwd=.9)

#### Legend

legend('topleft',legend = c(paste('Adj. r^2 ',substr(summary.lm(lm1)[9],1,4),sep=''),paste('Slope ',round(lm1$coefficients[2],2),sep=''),paste('Intercept ',round(lm1$coefficients[1],2),sep='')),text.col=c('lightseagreen','lightseagreen','lightseagreen'), bty='n')

#### Shut device

dev.off()

# Write out model parameters and total deviance as a .csv file

params = data.frame(learning.rate = lr, tree.complexity = tc, train.fraction = tf, n.trees = nt, param.set = sprintf('%02i',i), op.tree.num = brt.gbm.step$n.trees, test.deviance = total.deviance, train.deviance =  brt.gbm.step$cv.statistics$deviance.mean)

write.csv(params, file=paste('ul.maxtemp.model.summary_',sprintf('%02i',i),'.csv',sep=''),row.names=F)

### Make some plots from training data
#### Start with obs versus preds

### Start .png driver

png(paste(output.dir,'Training_Data_Obs_Versus_Preds.png',sep=''),units='cm',height=12, width=12, res=1000)

#### Plotting limits

lims = range(c(train.data[,3],train.data[,21]),na.rm=T)

#### Plot up testing obs versus preds

plot(train.data[,3],train.data[,21], xlab = 'Obs Training Data', ylab = 'Training Fitted Values', main = paste('Obs vs Preds - Deviance - ',round(brt.gbm.step$cv.statistics$deviance.mean,4),sep=''), xlim=lims, ylim=lims, col='olivedrab1')

#### 1:1 line

abline(a=0,b=1, lty=2, lwd=1, col='red')

### LM

lm1 = lm(train.data[,21]~train.data[,3]) # Perform a linear model using the Empirical Data as y-values and the preds (AWAP or microCLIM) as x-values

#### Plot lm relationship

abline(lm1,col='lightseagreen',lty=2, lwd=1)

#### Legend

legend('topleft',legend = c(paste('Adj. r^2 ',substr(summary.lm(lm1)[9],1,4),sep=''),paste('Slope ',round(lm1$coefficients[2],2),sep=''),paste('Intercept ',round(lm1$coefficients[1],2),sep='')),text.col=c('lightseagreen','lightseagreen','lightseagreen'), bty='n')

#### Shut device

dev.off()

#### Now plot obs versus resid

### Start .png driver

png(paste(output.dir,'Training_Data_Obs_Versus_Resids.png',sep=''),units='cm',height=12, width=12, res=1000)

#### Plotting limits

xlims = range(train.data[,3],na.rm=T)
ylims = range(train.data[,3] - train.data[,21],na.rm=T)

#### Plot up testing obs versus preds

plot(train.data[,3],train.data[,3]-train.data[,21], xlab = 'Obs Training Data', ylab = 'Training Resid Values', main = paste('Obs vs Preds - Deviance - ',round(brt.gbm.step$cv.statistics$deviance.mean,4),sep=''), xlim=xlims, ylim=ylims, col='olivedrab1')

### LM

lm1 = lm(train.data[,3]-train.data[,21]~train.data[,3]) # Perform a linear model using the Empirical Data as y-values and the preds (AWAP or microCLIM) as x-values

#### Plot lm relationship

abline(lm1,col='lightseagreen',lty=2, lwd=1)

#### Legend

legend('topleft',legend = c(paste('Adj. r^2 ',substr(summary.lm(lm1)[9],1,4),sep=''),paste('Slope ',round(lm1$coefficients[2],2),sep=''),paste('Intercept ',round(lm1$coefficients[1],2),sep='')),text.col=c('lightseagreen','lightseagreen','lightseagreen'), bty='n')

#### Shut device

dev.off()

### Write out test data

test_out = data.frame(test.data,ulmaxpreds=gbm.predicted)

write.csv(test_out,file=paste(output.dir,'UL_BRT_Test_Data_Plus_Preds.csv',sep=''),row.names=F)

### Done